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Abstract 

We introduce a mixed discontinuous/continuous finite element pair for ocean mod- 
elling, with continuous quadratic layer thickness and discontinuous velocity. We in- 
vestigate the finite element pair applied to the linear shallow-water equations on an 
/-plane. The element pair has the property that all geostrophically balanced states 
which strongly satisfy the boundary conditions have discrete divergence equal to 
exactly zero and hence are exactly steady states of the discretised equations. This 
means that the finite element pair has excellent geostrophic balance properties. We 
also show that the element pair applied to the non-rotating linear shallow-water 
equations does not have any spurious small eigenvalues. We illustrate these proper- 
ties using numerical tests and provide convergence calculations which show that the 
numerical solutions have errors which decay quadratically with element edge length 
for both velocity and layer thickness. 



1 Introduction 



A number of finite element pairs have been proposed for the rotating shallow- 
water equations, including the PInc ~ PI and PI — iso P2 — PI elements 



(investigated and compared to several other element pairs in (Le Roux et al 



1998)), the RTO elements (introduced in Raviart and Thomas (1977) and 
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proposed for the shallow- water equations in Walters and CasuUi (1998)) and 
equal-order elements with stablisation (also proposed in Walters and CasuUi 



(1998)); all of these elements have been shown to perform well when inte- 
grating the rotating shallow-water equations. In this paper we investigate the 
numerical properties of the P1dg-P2 finite element pair applied to the linear 
shallow-water equations on an /-plane in order to investigate the suitability 
of the element for shallow-water ocean modelling. The finite element pair con- 
sists of discontinuous linear elements for velocity and continuous quadratic 
elements for layer thickness. Even though the layer thickness has shape func- 
tions which are one order higher than velocity, there are still more degrees 
of freedom in the space of discontinuous linear functions than the space of 
continuous quadratic functions (except in meshes with very few elements), 
which is a necessary (but not sufficient) condition for the absence of spurious 
pressure modes (modes with high spatial frequency but small eigenvalues in 
the wave operator which can pollute the numerical solution with noise). In 



Cotter et al. (2008), it was shown that in one dimension, this choice leads to a 
discretisation of the wave equation without rotation which does not have any 
spurious modes. It was also shown that the dispersion relation is very accurate 
for the first half of the discrete spectrum. A number of numerical tests were 
carried out on two- and three-dimensional unstructured meshes which showed 
that there were no spurious eigenvalues present. In fact, as we show in this 
paper, the discretised Laplacian obtained by combining the first-order discrete 
divergence and gradient operators is the same as the usual Galerkin finite el- 
ement discretised Laplacian obtained by multiplying the Laplacian by a test 
function and integrating by parts. Consequentially, the discretised equations 
without rotation do not have any spurious eigenvalues, also the element pair 
applied to the incompressible Navier-Stokes equations leads to an LBB-stable 
discretisation without spurious pressure modes. For a general discussion of 



LBB stability, see Gresho and Sani 


(2000 


); 


Karniadakis and Sherwin 


(2005) 


Karniadakis and Sherwin 


(2005) also contains an exposition of the discontinu 



ous Galerkin method. For applications of the discontinuous Galerkin method 



to waves equations see Ainsworth et al. (2006), and for some applications of 



the discontinuous Galerkin method to the rotating shallow-water equations 
see lAmbati and Bokhovel (|2007[); iLevin et all (|2006l); iBernard et all (|2007|); 



Giraldo (2006). 



In this paper we concentrate on the interaction of the geostrophic modes 
with the inertia-gravity waves which is crucial to the good representation of 
large-scale dynamics. We find that not only does the finite element pair allow 
for accurate representation of geostrophically-balanced states, these states are 
completely uncoupled from the inertia-gravity waves: the states are exactly 
steady as in the unapproximated partial-differential equations. In section 2 
we introduce the element pair applied to the linear shallow water equations, 
show that the element pair has a discrete Laplacian without spurious eigen- 
values, and show that the element pair has exactly steady geostrophic states. 
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In section 3 we verify these results with numerical tests. We also show nu- 
merical calculations using Kelvin waves which arc gcostrophically balanced 
in one direction; these waves are a good test of preservation of balance. The 
results do not show any radiating inertia-gravity waves. We provide conver- 
gence test results using the Kelvin wave exact solution which confirm that 
the errors spatial discretisation converges quadratically for both velocity and 
layer thickness, indicating that the element pair is stable. Finally we give a 
summary and outlook in section 4. 



2 The mixed element 



In this section we describe our mixed element formulation applied to the linear 
shallow-water equations on an /-plane. 



2.1 Mixed continuous/discontinuous Galerkin discretisation 

We start with the linearised shallow-water equation on an /-plane in non- 
dimensional units 

Ut+ r^k X U+ :^Vh = 0, U= {Ui,. . . ,Ud), (1) 
Ko i'r 

/it + V • u = 0, (2) 
where u is the velocity, h is the perturbation layer thickness, k is the unit 



vector in the 2;-direction, Ro = U / fL is the Rossby number, Fr = ^^U /gH is 
the Froude number, ?7 is a velocity scale, L is a horizontal length scale, H is 
the mean layer thickness, / is the Coriolis parameter and g is the acceleration 
due to gravity. The boundary conditions are 

u - n — Q on do, (3) 

where dVt denotes the boundary of the domain fi, and n is the normal to 
do.. To obtain the discontinuous/continuous Galerkin form of the equations 
we multiply equation (1) by a discontinuous test function w and equation (2) 
by a continuous test function (j) and integrate over an element E to obtain 



d 

di 



[ w -udV + f w ■ k X udV^ — ^[w-VhdV. (4) 
Je Ro Je Fr"^ Je 

f (f)hdV^- I (t)V -udV. (5) 
J e Je 



d_ 
dt 
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Wc then integrate equation (5) by parts, and make use of the boundary con- 
ditions (3) to obtain 



/ w -udV + I w - kxudV ^ — \y I w-VhdV. (6) 

JE Ro Je ¥r Je 

f (l)hdV= I V(t)-udV (7) 

Je Je 

- I n-u(i)dS, (8) 

JdE\dn 



d 

dt 



where ii is the value of u on the element boundary dE, determined by the 
particular choice of discontinuous Galerkin scheme which is chosen (the value 
on the upwind face, for example), and where n is the outward-pointing unit 
normal to the surface dE. Conservation requires that u takes the same value 
on either side of each face. We sum these equations over all elements and the 
surface terms cancel since is continuous. This gives the form of the equations 
that we will discretise: 



di 



f w - udV + f w ■ k X udV = — \ [ w - VhdV, (9) 

Jn Ro Jn Fr^ Jn ^ ^ 

[ (j)hdV^ [ V(l)-udV. (10) 
Jn Jn 



d 

di 



Derivatives are only applied to the scalar functions h and (f) and not the vector 
functions u and w which we shall discretise with discontinuous elements. To 
add nonlinear advection it is necessary to develop surface integrals on the 
boundaries of the elements following the standard discontinuous Galerkin finite 
element approach. 



2.2 The P1dg-P2 element 



In this subsection we develop the P1dg-P2 discretisation for the shallow-water 
equations. We make the choice that u and w are approximated by discontinu- 
ous linear finite element functions and lu'', whilst 4> and h are approximated 
by continuous quadratic linear finite element functions and (j)^ . 

The Galerkin finite element approximation of equations (9,10) is then 

^ I ■u^dV+^ I -kx u^dV = — ^ / -Vh^dV, 
dt Jn Ro Jn Fr Jn 

^ I ^^h^dV= I V4)^-u^dV, 

dt Jn Jn 
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u node Q 



h node 



Fig. 1. Figure showing the distribution of nodes in the two-dimensional P1dg-P2 
element. Each element contains three nodes for each of the two components u and 
V of velocity, and six nodes for the layer thickness, but the latter nodes are shared 
across element boundaries since the layer thickness space is continuous. 

for all test functions (p^ and in the specified spaces. 



2.3 Properties of discretised Laplacian 

The first property to note for the P1dg-P2 element pair is that the discretised 
gradient in the PIdg velocity space of a P2 function obtained from 



at each point. To prove this, note that all continuous P2 functions h have 
PIdg gradients. This means that we may choose = — V/i^ and hence 



Since and Vh^ are piecewise polynomials this means that they are identi- 
cally equal. 

The discretised Laplacian in the layer thickness space is obtained by ap- 
plying the discretised divergence to the discretised gradient q^: 




for all PIdg test functions w^, satisfies 




5 



Jci Jn Jdii 

Jn Jan on 



for any P2 test function 0"^, which is the standard Galerkin finite element dis- 
cretisation of the Laplace operator obtained by multiplying by a test function 
and integrating by parts. The properties of this operator are well-known in 
the finite element literature; in particular it has no spurious eigenvalues. 



2.4 Exactly steady geostrophic modes 

In the linear shallow-water equations with rotation, the geostrophic balanced 
modes with 

Ut = ^ u = V^ip, (11) 

are steady in time, where 

Ro 

This is because ht = since V ■ tt = for these modes. In this section we show 
that the balanced states in discretisations with the P1dg-P2 element pair are 
also completely steady; this means that the P1dg-P2 element pair represents 
balanced states very well and so is ideal for shallow-water ocean modelling. 

The geostrophically balanced states in the finite element discretisation satisfy 

[w^-u^dV^O, =^ [ -u^dV ^ [ -V^i/j^dV, (12) 
Jn Jn Jn 

for all Pl^*^ test functions w^. 



d 

di 



The property described in the previous section can be trivially extended to 
show that the finite element velocity obtained from this equation for a given 
■0'' satisfies equation (11) with u — and ■0 = V'*- We can use this property 
to show that any geostrophically balanced velocity field u obtained from a 
streamfunction ip which is constant on the boundary satisfies the discrete 
divergence equation 

- / W(t)^ -u^dV = 0, 
Jn 

for all P2 test functions (f)^ . To prove this, note that 
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- [ V(f)^ ■u^dV=- [ -V^ip^dV (13) 
= Y. I (t^^^ -^^^^ -Y. I 0^n-VVd5 (14) 

= -y I (t)^\\n-V^^l)^]]dV - I -71(15 = 0(15) 

P Jr ^ ' Jan ""^^ 

=0 -0 

where Y1,e indicates a sum over all elements E, dE is the boundary of element 
E, X^r indicates a sum over all orientated internal element boundaries in the 
mesh, and [[/]] indicates the jump in a function / across a surface V. In 
equation (14) the normal component of velocity vanishes exactly on dVL as h 
is constant on dVt and the balanced velocity is obtained from the pointwise curl 
of the streamfunction ip. In (15) the jump in the normal component of V'^ip 
vanishes because the tangential derivative of functions in P2 is continuous 
across element boundariesPH 

The proof of this property is easily extended to the general PnDG-P(n+l) 
element pair i.e., rath order discontinuous velocity and (n+ l)-th order contin- 
uous layer thickness. It is also easily extended to the three-dimensional case 
in which u = V A for any vector field ^ which is constant on the boundary. 



3 Numerical tests 



In this section we illustrate and explore the properties of the P1dg-P2 element 
applied to the linear rotating shallow-water equations. 



3.1 Representation of geostrophic balance 



Le Roux et al. (1998) tested a number of element pairs for their ability to 
represent geostrophic balance. This was done by selecting a streamfunction 
field, computing the balanced velocity field from equation (12), and plotting 
streamlines. Element pairs were compared by the smoothness of the resulting 
streamlines on structured and unstructured meshes. Here we just note that, 
as described in the previous section, the balanced velocity for the P1dg-P2 
element is obtained from the pointwise gradients of streamfunction and so 
streamlines of the discretised balanced velocity field are simply contours of 



^ The right-hand side of (13) can be shown to vanish for general functions from 
the space (which contains the P2 functions) by taking a convergent sequence of 
smooth functions and passing to the hmit. However, the extra property of continuous 
tangential derivatives of Pn functions facilitates the simpler proof given here. 
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Fig. 2. Left: Streamlines of balanced velocity obtained from a Gaussian stream- 
function distribution using the P1dg-P2 discretisation. The streamlines are very 
smooth showing that the discretisation does not introduce spurious oscillations. 
This is because in this case the balanced velocity can be obtained by taking the 
pointwise (strong) gradients of the streamfunction. Right: The mesh used for this 
calculation. The mesh was deliberately distorted to illustrate that this property is 
not dependent on mesh quality. 



the discretised streamfunction field. This means that the balanced velocity 
field is actually as accurate as possible for the P2 streamfunction field. Plots 
of some resulting streamlines are given in figure 2; for comparison with other 
element pairs see Le Roux et al. (1998). 



3.2 Steady states 



In Le Roux et al. (1998), another numerical test was performed in which 



the linear rotating shallow-water equations were initialised in a geostrophic 
state; streamlines were plotted after some time which showed that the PI iso 
P2— PO — 3 element pair (proposed in that paper) preserved the steady state to 
excellent accuracy. In the case of the P1dg-P2 element, we have already shown 
in the previous section that geostrophic states are exactly steady so it remains 
to verify this numerically. Using the mesh shown in figure 2, we computed ran- 
domly generated streamfunction fields with ijj = on the boundary together 
with their geostrophically balanced velocity fields obtained from equation (11), 
and integrated the equations in time using the Crank-Nicholson method. We 
observed that the layer thickness h and velocity u remained constant up to 
machine precision, confirming that the geostrophic modes are completely un- 
coupled from the inertia-gravity waves. For the time evolution of geostrophic 
states using other element pairs, see Le Roux et al. (1998). 
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Fig. 3. Plot showing the mesh used for Kelvin wave tests. 
3.3 Kelvin waves 

We tested the P1dg-P2 element using a Kelvin wave initial condition; the 
Kelvin wave is a trapped coastal wave which is geostrophically balanced in 
the direction normal to the coast which propagates at the fast gravity wave 
speed 1 /Fr for the case of a straight coastline. The aim of the test is to verify 
that the Kelvin wave does not shed any spurious inertia-gravity waves. We 
used the circular Kelvin wave initial condition given by 

/i(r, e) = e('-''°)/^°cos^, 
Fr 

Ur — 0, 

with Ro = 0.1 and Fr = 1. The Kelvin wave propagates around the circular 
coast, maintaining geostrophic balance in the normal direction. The mesh 
used for the discretisation is shown in figure 3. We integrated the equations 
in time for > i > 100 using the Crank-Nicholson method and a time step 
size At = 0.01. Figure 4 shows the layer thickness at various times: there are 
no spurious gravity waves observed, which means that the P1dg-P2 element 
pair is maintaining geostrophic balance in the normal direction as well as the 
Kelvin wave structure. 

To check convergence of the method we integrated a Kelvin wave in the rect- 
angular domain = {a;:— 15<x<15,0<|/<3} with initial condition 

h = e-^/^°e-(=^-^)' , u = (e-2//i^°e-(^-^)' , 0) . 
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Fig. 4. Plots showing contours of layer thickness h at times t =0 (top left), 30 

(top right), 60 (bottom left) and 90 (bottom right) for the circular Kelvin wave test 
case. No spurious oscillations are observed, which verifies that the P1dg-P2 element 
maintains geostrophic balance in the normal direction. 



If this initial condition is used in the domain fl°° = {cc:— cx3<a;<cxD,0< 
y < oo}, then the equation has the exact solution 

h = e-^/^°e-("+*/^'^'-^)', u = (e-^/^°e-(^+*/^'-^^', 0). 

We integrated the system to time t — 10. For this time interval the solution 
is almost zero for y > 1 and |x| > 6 and so the exact solution is a good 
approximation. The timestep was chosen to have a wave Courant number 
of less than 0.1 for all simulations so that the errors are dominated by the 
spatial discretisation. We refined the mesh isotropically in space in the region 
where the solution was non-zero during the calculation and computed the 
L2 errors of the velocity and the layer thickness for various element edge 
lengths in the refined region. Plots of the numerical errors are given in figure 
5. A linear regression on these values showed that the velocity errors were 
proportional to 2.19 and the layer thickness errors were proportional to 1.98. 
These results suggest that the errors in velocity and layer thickness in the 
spatial discretisation scale quadratically with the edge length, as would be 
expected from approximation theory. They are also an indication that the 
element pair is stable: if the element pair were unstable then there would be 
spurious modes present which would lead to slower convergence than that 
expected from approximation theory. 
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Fig. 5. Convergence plots for tests with a Kelvin wave propagating along a flat coast 
performed on unstructured isotropic triangular meshes, showing error e against 
element edge length Ax. Left: L2 error in velocity plotted against element edge 
length. Right: L2 error in velocity plotted against element edge length. Both plots 
show that the errors scale with Ax'^ as Ax — > 0. 

4 Summary and Outlook 



In this paper we introduced the P1dg-P2 element pair applied to the linear 
shallow-water equations on an /-plane. We showed that the element pair has 
the property that all geostrophically balanced states which strongly satisfy 
the boundary conditions are exactly steady since their discrete divergence is 
identically zero. This means that the element pair has excellent geostrophic 
balance properties. We verified these properties by computing the evolution of 
balanced states, and by simulating Kelvin wave solutions which are geostroph- 
ically balanced in one direction. Finally we gave convergence test results which 
show that the numerical solutions have errors which decay quadratically with 
element edge length; this verifies the LBB-stability properties discussed in 



Cotter et al. (2008) 



In future work we shall compare this element pair with other low-order ele- 
ment pairs such as the POdg ~ -PI; P^nc ~ and RTO pairs. Whilst the 
discontinuous velocity means that the P1dg-P2 pair has a large number of 
degrees of freedom per element, the remarkable accuracy of the first half of 



the dispersion relation (noted in Cotter et al. , 2008 ) suggests that the element 



may be competitive, especially given its excellent treatment of geostrophic bal- 
ance, and local conservation of momentum. The higher-order extensions such 
as P2£,G ~ P3 will also be examined. We shall investigate the performance of 
the element once nonlinear advection has been introduced. 



A key advantage of this element pair is that the extension to three dimen- 
sions is also LBB-stable; the property that geostrophically balanced states are 
exactly divergence-free also extends to the three dimensional case. We shall 
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investigate the performance of this element pair in fully three-dimensional 
unstructured mesh ocean modelling in the ICOM model (Pain et al. , 2005). 
We also expect that if the buoyancy is discretised using PI^g elements, then 
the discretisation will also preserve hydrostatic balance very well; this will be 
investigated in future work. 
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